NASA Technical Memorandum 102800 


Extension of a Streamwise 
Upwind Algorithm to a 
Moving Grid System 

Shigeru Obayashi 

MCAT Institute, San Jose, California 

Peter M. Goorjian and Guru P. Guruswamy 
Ames Research Center, Moffett Field, California 


April 1990 


IWNSA 

National Aeronautics and 
Space Administration 

Ames Research Center 

Moffett Field, California 94035-1000 




EXTENSION OF A STREAMWISE UPWIND ALGORITHM 
TO A MOVING GRID SYSTEM 


Shigeru Obayashi,* Peter M. Gooijian, and Guru P. Guruswamy 
Ames Research Center 


SUMMARY 


A new streamwise up jvind algorithm has been derived to compute unsteady flow fields with the use 
of a moving-grid system. The temporally nonconservative LU-ADI (lower-upper-factored, alternating- 
direction-implicit) method has been applied for time-marching computations. A comparison of the tempo- 
rally nonconservative method with a time-conservative implicit upwind method indicates that the solutions 
are insensitive to the conservative properties of the implicit solvers when practical time-steps are used. Us- 
ing this new method, computations have been made for an oscillating wing at a transonic Mach number. 
The computed results confirm that the present upwind scheme captures the shock motion better than the 
central-difference scheme based on the Beam -Warming algorithm. The new upwind option of the code 
allows larger time-steps and thus is more efficient, even though it requires slightly more computational 
time per time-step than the central-difference option. 

INTRODUCTION 


A code, ENSAERO, is being developed at Ames using the Euler/Navier- Stokes equations for com- 
puting the unsteady aerodynamics and aeroelasticity of aircraft. The capability of the code has been demon- 
strated by computing vortical and transonic flows over flexible swept wings (refs. 1 and 2). The flow fields 
were calculated by a time-accurate, finite-difference scheme based on central differencing. 

The purpose of this study is to enhance the algorithm capability of the present code. In this re- 
spect, the use of a new upwind scheme in comparison to the current central-difference (CD) scheme is 
investigated. The CD scheme requires an artificial dissipation to stabilize computations. Such artificial- 
dissipation models lead to more dissipative solutions than upwind schemes. In addition, the CD scheme 
is sensitive to the amount of dissipation and needs a specific dissipation coefficient for each case. On the 
other hand, upwind schemes do not require that any coefficient be specified. 

Recently, a streamwise upwind algorithm has been developed and applied to treat steady-state prob- 
lems of transonic flows over wings (ref. 3) and vortical flows over a delta wing (ref. 4) on fixed grids. The 
main feature of the streamwise method is the use of the local stream direction, flow velocity, and pressure 
gradient. The switching of flux evaluations always takes place at sonic values, where shock waves may ex- 
ist. Therefore, this method fellows the flow physics more closely than conventional upwind methods based 
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on dimensional splitting. The computed results confirm the higher resolution of the present algorithm over 
the CD scheme, as well as over other upwind schemes. 

In this paper, the stream wise upwind algorithm has been extended from fixed coordinates to moving 
coordinates for computing flows over moving components. The streamwise upwind scheme was applied to 
steady-state problems by using the lower-upper-factored, altemating-direction-implicit (LU-ADI) method 
in order to accelerate convergence (ref. 4). The same LU-ADI method, which is first-order accurate but 
nonconservative in time, is used in the present unsteady computations. In order to check the validity of the 
LU-ADI method for computations over moving grids, a conservative implicit version of the streamwise 
upwind scheme is also considered. The resulting algorithm has been implemented in the code as a finite- 
volume, upwind option, in addition to the previous CD option. The updated code has been successfully 
applied for computing unsteady transonic flows over an oscillating wing. The computed unsteady pressures 
are compared with the experimental data. 


GOVERNING EQUATIONS 


The nondimensionalized thin-layer Navier-Stokes equations used in this study can be 
conservation-law form in a generalized body-conforming curvilinear coordinate system ( £, 
follows: 


drQ + d^E + dr,F + d<G = — d<;G v 

ilC 


written in 

V, C ) as 
( 1 ) 


The Euler equations are obtained by setting the viscous flux vector G v equal to zero. The vector of con- 
served quantities Q and the inviscid flux vectors E, F, and G are 
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The Cartesian velocity components u, v, and w are nondimensionalized by the free-stream speed of sound 
Coo", the density p is nondimensionalized by the free-stream dejnsity Poo! ar *d the total energy per unit volume 
e is nondimensionalized by p x>c^,. The viscous flux vector G v is given by 


with 
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where Re is the Reynolds number, Pr is the Prandtl number, c is the speed of sound, and J is the trans- 
formation Jacobian. Pressure is related to the conservative flow variables Q, through the equation of state 


for a perfect gas: 

P = (7 - 1) [ e - ^-(u 2 + v 2 + w 2 ) j 


( 2 ) 


where p is the fluid density a id e is total energy per unit of volume of the fluid. 

For the inviscid case, the viscous flux G v is replaced by 0. For the viscous case, the viscosity 
coefficient p in G v is computed as the sum of p* + p t where the laminar viscosity p* is taken from the 
free-stream laminar viscosity, assumed to be constant for transonic flows, and the turbulent viscosity p< is 
evaluated by the Baldwin-Lomax model. 5 


NUMERICAL ALGORITHM 


The space-discretized form of equation (1) can be written as 
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where a second-order central -difference evaluation is applied to the viscous term. 

The evaluation of the inviscid fluxes is based on the finite-volume cell-centered scheme. To be 
consistent with the finite-difference scheme in ENSAERO, the metrics are defined at each grid point where 
the flow variables are stored. The surface vector of cell-interface, which is necessary for the finite-volume 
formulation, can be obtained by averaging the metrics at the adjoining points. The free-stream preservation 
of this metric evaluation was shown in reference 6. 
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Streamwise Upwind Algorithm 


The formula of the streamwise upwind algorithm on the fixed grid system (& = Vt - Ct - 0 ) can 
be written in the cell-interface flux for the ^-direction with a surface vector S normal to the £-( plane as 
follows: 

F(Qi,Q r , S ;+ i) [1 + sign (Vi) cos 2 + s t A*(pg), e 3l cos 2 9i 
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and where q is the velocity magnitude, e» = ( 1 , u, v, w, H) T , t d - (0 , k x , k y , k z , k x u + k y v + k z w) T , 
A • = - r - -i, V = k x u + k y v + k z w, k x = r) x /\ Vtj| , and so on. The basic scheme is first-order accurate with 
l = j and r = j + 1 . The variables are averaged unless defined as a difference between left and right states. 
The switches si >r and the rotation angle 6 will be defined later. The cos 2 <? terms represent the projection on 
the coordinate axis from the streamwise flux vector splitting, and the sin 2 0 terms represent the projection 
on the crossflow plane for those terms that use the flux difference splitting. More details for the formulas 
are given in reference 4. 


To extend equation (4) from a fixed-grid system to a moving-grid system, the flow velocity relative 
to the fixed grid, q = (u, v,tu), is redefined as the flow velocity measured relative to the moving grid, 
q = (u — %t,v — yt,w — zt) ■ This is consistent with the modification of the definition of the contravariant 
velocity from a fixed-grid system, 

V fixed = rj x u + ri y v + rj z w 
= V 7] • q 


to a moving-grid system, 

Vmoving = rjt + rfcu + riyV + r) z w 

= rj x (u — x t ) + r\ y (v - y t ) + r] z (w — z t ) 

= Vtj • q 

where rj t = —T] x Xt — TjyPt ~~ Vz z t ■ Thus, equation (4) can be rewritten for a moving-grid system by using 
the following redefined variables: q 2 = (u — x t ) 2 + (v — y t ) 2 + (w — z t ) 2 and V - k x (u — x t ) + k y (v — 
yt) + k z ( w — z t ) . (Note that e s and e d remain unchanged because they originate in the eigenvectors of the 

Jacobian dF/dQ.) 
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Equation (4) is written in a vector form. If the formula is rewritten in component form, the present 
algorithm can be summarized with a surface vector S and a motion of its centroid, x t = ( x t ,yt, z t ) > as 
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where q and V are defined above, and where Ai = (c m — | Vm|) A 2 = {c m — |V r m|)PmAV, and 
kt = -k x x t - k v y t - k z z t . The averaged state (m) is defined for p, u, v, w, and H by the arithmetic 
average of the left (J) and right (r) states. 

The switches si and i r are defined in the manner of Godunov’s method as follows. For V > 0, 
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and M = q/c. 

A simple way to evaluate the rotation angle is to use cos0 = V /q. In supersonic flow fields, 
however, it is important to detect whether the velocity projected to the grid line is beyond the Mach cone. 
Thus, V/q is replaced by M ■ V /q = V/c. If V /c becomes larger than one, cos0 is set to one. To avoid 
expansion shocks, the rotation angle is determined by a mixture of averaged (m) and pointwise (l, r ) values. 
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The following relation is used here for evaluating $€[0,1] because of the smoothness: 


$ = max 



0 


( 8 ) 


where p\ and p 2 denote upstream and downstream pressures, respectively. The sine is determined by an 
arithmetic average of the cosines: sin 2 0 = 1 — j(cos 2 0j + cos 2 0 r ). 

Higher-order schemes are constructed from a one-parameter family, «, of interpolations of the 
primitive variables, p, u, v, w, and p. For example, 

P< = {l + ^[(l-«)V+(l + /0A]}p ; 

Pr = {! — 1 + k) V + ( 1 - «)A]}p ;+1 (9) 


where V and A are backward and forward difference operators, respectively (ref. 7). For the third-order 
scheme, k = j, Koren’s differentiable limiter (ref. 8) is used here. The limiter $ is calculated as 
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where a small constant e, e = 10 _6 typically, is added to prevent the division by zero. The same formulas 
are used for the other primitive variables. 


LU-ADI Method 


One of the time-marching methods for the present upwind scheme is the LU-ADI factorization 
method proposed by one of the present authors (ref. 6). The LU-ADI method is a compromise between 
ADI and LU factorization. Applied to equation (3), this method is written as, 

{T^aDaUaT^) x (T^LbDbUbT- 1 ) x ( T^LcDcUcT, f 1 ) x AQ n = AtR n (11) 

where R represents the right-hand side of equation (3). The change of volume in time is neglected, assum- 
ing translational movements of the grids. On the left-hand side, the original ADI operator of the Beam- 
Warming method is rewritten using the diagonal form (ref. 9) and the first-order accurate Steger- Warming 
flux-vector splitting (ref. 10). For example, in the p -direction, 

I + A t8fjB 

= T v ( I + AtV v D + B + At\D- B )T~ x 

= T n U~ AtDg\j + A tV n D + B HI + A t\D B \})~\l + AtD + B \j + AtA^DpT 1 - 1 

= T v L b DbUbT- x (12) 

This factorization is the approximate LDU (lower-diagonal-upper) factorization. This is more stable than 
simple LU factorization since the diagonal element always has the absolute value of the eigenvalues. 
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Approximate Block ADI Method 


The LU-ADI method described in the previous section is nonconservative in time owing to the 
diagonalization. To investigate the significance of this temporal nonconservativeness, an alternative ap- 
proach is considered here. A time-conservative method can be constructed using a block- tridiagonal solver 
similar to the Beam-Warming method. Since true Jacobians of the numerical fluxes of the present upwind 
algorithm are expensive to compute, approximate Jacobians are used here. 


To construct an implicit method for the present upwind algorithm, it is easier to start from the vector 
form (eq. (4)) rather than from the component form (eq. (5)). From equation (4), the first-order-accurate 
flux can be rewritten as 
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To simplify the formula further, the e 4;± i term in the second line of the above formula is replaced by e s; . 
Then, the Jacobian dF ± /dQ can be approximated as 
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and where V f = k x u + k y v + k z w. The switches and the angles are evaluated identically to the explicit 
formulas. Similar to the results discussed in reference 1 1 , this approach yields a method that is less stable 
than the LU-ADI method. 


Neglecting the viscous terms, the approximate block ADI method can be written as 
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In contrast to the LU-ADI method, the resulting method is conservative. However, both methods are first- 
order accurate in time. This method requires twice as much CPU time as the LU-ADI method. 


RESULTS 


Time-accuracy is an essential requirement for aeroelastic computations. Numerical schemes used 
for flow calculations in aeroelasticity must guarantee the correct calculation of the amplitude and phase 
of unsteady pressures. In order to verify the time-accuracy of the present code, unsteady flows over a 
rectangular wing undergoing a prescribed oscillatory motion are computed. 

The test cases consider inviscid and viscous unsteady flows over a rectangular wing with a NACA 
64A010 airfoil section and an aspect ratio of 4. Steady and unsteady measured data from wind-tunnel tests 
of this wing are given in reference 12. The unsteady data are given for the case in which the rigid wing 
is oscillating in the pitching mode, a(t) = ot m — aosin(u>t), about an axis at x/c r = 0.5, where Cr is 
the reference chord length, and u> is the pitching frequency in radians per second. The inviscid flow is 
computed at M m = 0.8 with a mean angle of attack a m = 0°, pitch amplitude a 0 = 1°, and reduced 
frequency k= 0 .27 (k= uj c r j Uoo ) • la addition, for the viscous flow, the Reynolds number based on the 
wing chord is Re = 2 x 10 6 . The Baldwin-Lomax model (ref. 5) is used to compute the turbulent eddy- 
viscosity coefficient for this case. 

The C-H-type grid system is used for the present computations. The £, rj, and C coordinates repre- 
sent the chordwise, spanwise, and normal (to the wing surface) directions, respectively. For the inviscid 
computations, coarse and fine grids are used. The coarse grid contains 91 x 25 x 25 points in the £, rj, and 
^-coordinate directions, respectively, and the fine gnd contains 151 x 25 x 34 points. These grids have 
a minimum spacing normal to the wing of 0 .01 c r . The viscous grid also contains 151 x 25 x 34 points 
and has the minimum spacing normal to the wing of 0 .00005 c r . This yields a value of y + <6 at the first 
grid point above the wing surface. The grids are generated algebraically. 

Unsteady computations are started from the corresponding steady-state solution. The stability of 
unsteady computations is verified by increasing the numbers of time-steps per cycle until no change is 
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noticed in the unsteady results The convergence of the unsteady computations to a periodic flow is verified 
by comparing the results between cycles. For all cases considered here, the results for the third cycle give 
pressure profiles that are identical to those of the second-cycle results. Thus, the numerical transient is 
confirmed to disappear within two cycles. 

The present upwind results are also compared with the results obtained by using the existing CD 
method. Reference 13 discusses the CD results obtained from the same test case. The CD method uses 
the diagonal inversions of the Beam-Warming method (ref. 9), and thus it is first-order accurate but non- 
conservative in time. Its artificial dissipation model consists of second- and fourth-order dissipation terms 
controlled by the amount of the second-difference of pressure. 


Inviscid Solutions 


First, the inviscid coarse-grid (91 x 25 x 25 points) solutions are presented. Figure 1 compares the 
computed and measured unsteady, upper-surface pressure coefficients for the real and imaginary parts of 
the first Fourier component. Pressure coefficients are shown for various span locations. The computations 
were done by using the upwind algorithm with the LU-ADI method (UP-LU). The solid, dashed, and dotted 
lines in the figure indicate the results obtained using 1440, 720, and 360 time-steps per cycle of oscillation 
(steps/cycle), respectively. (The time-steps per cycle rate of 1440 corresponds to At ~ 0 .02 .) The solution 
profiles obtained using 1800 steps/cycle coincided with those using 1440 steps/cycle. Thus, the unsteady 
pressure profiles converged with respect to time-step sizes at 1440 steps/cycle. In addition, because the 
differences between the results with 720 and 1440 steps/cycle are not critical when making comparisons 
with experimental data, the computation with 720 steps/cycle is acceptable for numerical efficiency. 

Figure 2 shows the analogous CD results. In this case, the second- and fourth-order dissipation 
coefficients are fixed at 0.25 and 0.01, respectively (denoted as CD(0.01)). These coefficients are the 
values that were used in references 1, 2, 9 and 13. The CD results show less dependence on time-step 
size than do the UP-LU resulis, although the profiles are smeared out at the large-gradient region near the 
shock wave and leading-edge regions. In order to check the dependence of the CD method on the amount 
of numerical dissipation, computations were tried with half the amount of dissipation coefficients (0. 125, 
0.005). At 360 steps/cycle, the computation became unstable with these coefficients. At 720 and 1440 
steps/cycle, the computations were stable, but the resulting unsteady pressure profiles were still smeared 
out. 


Figure 3 shows the results using the approximate block ADI method applied to the present upwind 
algorithm (UP-BL). The UP-BL method was not stable for computations using 360 and 720 steps/cycle, and 
thus the computations were done using 1080, 2160, and 3600 steps/cycle. The solution profiles obtained 
with 3600 steps/cycle finally s howed good agreement with the UP-LU result converged at 1440 steps/cycle. 

Figure 4 compares the three methods, UP-LU, UP-BL, and CD(0.01), fixing the number of time- 
steps per cycle at 1440. Both UP-LU and UP-BL results give similar profiles at the shock motion, although 
the UP-LU method is nonconservative in time. Both upwind methods give better agreement with the 
experimental data in the region of the shock motion than the CD method. 
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These results indicate that the temporally nonconservative LU-ADI method is valid for unsteady 
computations, even with moving shock waves, when the number of time-steps per cycle is large (i.e., small 
At). In the present case, any number of time-steps per cycle greater than 720 will give a reasonable result. 
The CD method also uses the temporally nonconservative diagonal form. However, the CD result differs 
from both the upwind results. This indicates that the solution depends on the numerical dissipation more 
than on the time-conservative properties of the methods. 

The coarse-grid computations took 1 8.0, 19.5 and 42.3 /xsec per grid point per time-step for the CD, 
UP-LU, and UP-BL methods, respectively, on a CRAY-YMP computer using a single processor. Although 
the UP-BL method is the most expensive, its accuracy does not appear to compensate for the computational 
inefficiency. Thus, the UP-BL method will be dropped in the following computations. 

Next, the inviscid computations were done on a finer (151 x 25 x 34 points) grid. Since detailed 
comparisons with respect to time-step sizes are done on the coarse grid, the results are highlighted in the 
comparison of the UP-LU and CD solutions, as shown in figure 5. Both computations use 720 steps/cycle, 
as suggested in the coarse-grid case, because the time-step requirement for the fine-grid case is expected 
to be the same as that for the coarse-grid case (note that both grids have the same minimum spacing size 
normal to the wing). The fine-grid computations took 15.5 and 16.8 /rsec per grid point per time-step for 
the CD and UP-LU methods, respectively. 

For the CD method shown in figure 5, the two sets of the second- and fourth-order dissipation 
coefficients (0.25, 0.01) and (0.125, 0.005), indicated as CD(0.01) and CD(0.005), respectively, are used to 
illustrate the effect of the artificial dissipation. The CD solutions approach the UP-LU solution at 50% and 
77% semispan locations as the dissipation coefficients are reduced. However, the shock profile becomes 
too steep and the post-shock oscillation in the plot becomes deep at the root section. Further reduction of 
the dissipation coefficients (to 0.0625, 0.0025) made the computation unstable. Compared with the coarse- 
grid solution shown in figure 4, the CD solutions also tend to converge to the UP-LU solutions because of 
the grid refinement. Therefore, the UP-LU method is confirmed to give a less dissipative and thus more 
accurate solution than the CD method. 


Viscous Solutions 


To demonstrate the capability of the UP-LU method for viscous computations, the same case has 
been computed with the viscous terms. The viscous grid has the same number of grid points as the inviscid 
fine grid, but with a smaller spacing in the ^-direction normal to the wing. Because of the stiffness caused 
by the clustered viscous grids, time-step sizes used in the CD computations are usually determined by 
stability considerations rather than by time-accuracy (ref. 13). Two time-steps/cycle rates were chosen to 
demonstrate the robustness and accuracy of the UP-LU method. The results are again highlighted in the 
comparison of the UP-LU and CD solutions at specified time-step sizes. 

Figure 6 compares the viscous solutions using the UP-LU, CD(0.01), and CD(0.02) methods with 
720 steps/cycle. This number of time-steps per cycle is too small to obtain accurate and stable results for 
the CD method, a finding also reported in reference 13. The numerical instabilities in the CD results are 
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shown by oscillations in the pressure profiles of CD(0.01). However, the UP-LU result shows that the UP- 
LU method relaxes this stability requirement. Thus, time-step sizes for the UP-LU method are not limited 
by stability considerations even for the viscous case. 

Figure 7 shows the corresponding results with 1440 steps/cycle. This is the time-step rate sug- 
gested in reference 13. Here, the CD(0.01) solution becomes stable and agrees reasonably well with the 
experimental data. To check whether the amount of dissipation is reasonable, the dissipation coefficients 
were reduced by half, as was done for the inviscid case, but the CD(0.005) computation diverged for the 
viscous case. In figure 7, and analogous to the inviscid results shown in figure 5, the CD solutions approach 
the UP-LU solution as the dissipation coefficients are reduced. Again the UP-LU method is confirmed to 
give a less-dissipative and thus more accurate solution than the CD method. The viscous solutions in fig- 
ure 7 are similar to the inviscid solutions in shown in figure 5 because the experimental flow field does not 
contain strong viscous effect (ref. 12). 

These results also indie ate that the CD method is sensitive to both time-step size and to the amount 
of artificial dissipation added. Therefore, users of the CD method are required to find an adequate combina- 
tion of time-step size and dissipation coefficients on a case-by-case basis. In contrast, the UP-LU method 
is robust, and users do not have to find any dissipation coefficient. 

For the viscous case, the UP-LU computation requires 18.9 /isec per grid point per time step at 
a speed of 141 MFLOPS on a CRAY-YMP computer using a single processor, and the CD computation 
requires 17. 1 /isec at a speed of 132 MFLOPS. There is 1 1% increase of CPU time when the upwind option 
of the code is used. 


CONCLUSIONS 


A new stream wise upwind algorithm has been derived to compute unsteady flow fields with a 
moving grid system. The temporally nonconservative lower-upper-factored, altemating-direction-implicit 
(LU-ADI) method has been applied to compute flow over an oscillating wing at a transonic Mach number. 
A comparison of the temporally nonconservative method with a time-conservative version of the upwind 
scheme indicates that the solutions are insensitive to the time-conservativeness of the implicit solvers when 
practical time-step sizes are used. The temporally nonconservative upwind method was found to be more 
stable and twice as efficient as the time-conservative block ADI scheme, even though the block ADI scheme 
employed an approximate form of the Jacobians. Comparisons have been made between results obtained 
with both upwind schemes, with experimental measurements, and with computed results obtained using 
the existing central-difference method. 

Comparisons with experimental data show that the upwind algorithm predicts the shock motion 
better than the central-difference (CD) method. The CD solutions are also found to be sensitive to the 
amount of numerical dissipation, and as a result, the dissipation coefficients must be specified case by case. 
In comparison, the present upwind method does not require a dissipation coefficient. Thus, the method that 
combines the streamwise upwind and LU-ADI methods is proposed for practical computations. This new 
upwind method is incorporated in the aeroelastic code ENSAERO. The new upwind option of the code 
allows larger time-steps and thus is more efficient, even though it requires slightly more computational 
time per time-step than the CD option. 
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Figure 5. Comparison of in viscid upper-surface unsteady pressures computed by using the upwind LU- 
ADI method and the central-difference method with two sets of dissipation coefficients: 151 x 25 x 34 
grid. Moo = 0.8, <v ro = 0°,k=: 0.27, and a 0 = 1 °. 
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Figure 7. Comparison of viscous upper-; 
steps/cycle: 151 x 25 x 34 jpid, Re = 2 x 


face unsteady pressures computed by using 1440 time 


) 6 , Moo = 0.8, a m = 0°,k= 0.27, and a 0 = 1°. 
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